Построение солитонов КдФ при помощи суперпозиции преобразований Бэклунда

В.Э. Адлер. Классические интегрируемые системы, весенний семестр 2020

Построение n-солитонного решения уравнения КдФ

u=2(∂_x^2) log W((Subscript[y, 1]),...,(Subscript[y, n])),    (Subscript[y, j])=(e^(Subscript[X, j]))+(((-1))^(n-j))(e^(-(Subscript[X, j]))),    (Subscript[X, j])=(Subscript[k, j])x+4(Subscript[k, j]^3)t+(Subscript[δ, j]),    0<(Subscript[k, n])<...<(Subscript[k, 1]),

по схеме, описанной в 10_Backlund transformations.pdf .

(* число солитонов *)
n=10;

(* задаем случайные собственные значения и сортируем их по убыванию,
   также задаем случайные сдвиги *)
K=Sort[RandomReal[{1/2,3/2},n],#1>#2&]
d=RandomReal[{-5,5},n]

(* фазы *)
X[j_]:=K[[j]]x+4((K[[j]])^3)t+d[[j]]

(* затравочные решения *)
q[]=0;
Do[q[j]=-K[[j]]Tanh[X[j]],{j,n,1,-2}]
Do[q[j]=-K[[j]]Coth[X[j]],{j,n-1,1,-2}]

(* формула суперпозиции *)
q[i___,j_,k_]:=q[i]+(((K[[j]])^2)-((K[[k]])^2))/(q[i,j]-q[i,k])

(* вычисление решения *)
u=-2D[q@@Range[n],x];

(* компилированная функция, это ускоряет счет *)
cu=Compile[{{x,_Real},{t,_Real}},Evaluate[u]];

(* построение графика *)
cgru[t_,xmin_,xmax_]:=Plot[cu[x,t],{x,xmin,xmax},
PlotRange->{0,2.1((K[[1]])^2)},
PlotStyle->Blue,
Filling->0,
PlotPoints->50,
PlotLabel->"t= "<>ToString[t]];

(* анимация по t *)
Manipulate[cgru[t,-50,50],
{{t,0},-10,10,Appearance->"Labelled"}]

По такой схеме можно построить без особого труда 10–12 солитонов, (с вронскианами – не больше 4). При большем числе и здесь счет становится медленным, а также накапливаются ошибки в арифметических операциях. Если немного постараться, то можно дойти до 15, но дальше Mathematica справляется плохо. На курсе 2018, Д. Чукланов довел число солитонов до 40, с использованием Си.